Libraries & functions

First, I’m loading the required libraries & functions

suppressPackageStartupMessages(library(Seurat, lib.loc = "/software/Seuratv4/lib/")) # Seurat v4.4.0, library for single-cell analysis
suppressPackageStartupMessages(library(SeuratObject, lib.loc = "/software/Seuratv4/lib/")) # For compatibility purpose
suppressPackageStartupMessages(library(data.table)) # For writing DE gene file
suppressPackageStartupMessages(library(ggplot2)) # For plotting
suppressPackageStartupMessages(library(ggpubr)) # For grid plotting
suppressPackageStartupMessages(library(plotly)) # For interactive plots
suppressPackageStartupMessages(library(crayon)) # Just for bolding the console output :D

cat(bold("Seurat"), "version", as.character(packageVersion("Seurat")), "\n")
## Seurat version 4.4.0
cat(bold("SeuratObject"), "version", as.character(packageVersion("SeuratObject")), "\n")
## SeuratObject version 4.1.4
cat(bold("data.table"), "version", as.character(packageVersion("data.table")), "\n")
## data.table version 1.16.4
cat(bold("ggplot2"), "version", as.character(packageVersion("ggplot2")), "\n")
## ggplot2 version 3.5.1
cat(bold("ggpubr"), "version", as.character(packageVersion("ggpubr")), "\n")
## ggpubr version 0.6.0
cat(bold("plotly"), "version", as.character(packageVersion("plotly")), "\n")
## plotly version 4.10.4
# Color-blind friendly palette 3-colors
cbPalette <- c("#E69F00", "#000000", "#56B4E9")
steinPalette <- c("#80C980", "#BDADD4", "#376CB0", "#FBBF85", "#F0027E")
tissuePalette <- list(`Pan_neuro_control`="#E69F00", `Pan_neuro_ND75KD`="#56B4E9")
celltypePalette <- list(`unannotated`="#D3D3D3", `epithelial cell`="#E0AC69", `skeletal muscle of head`="#895129", `cone cell`="#FFFFB7", `outer photoreceptor cell`="#FFF060", `photoreceptor cell`="#F2DB00", `photoreceptor-like cell`="#F0C200", `adult optic chiasma glial cell`="#d896ff", `adult brain cell body glial cell`="#efbbff", `adult brain perineurial glial cell`="#BCBDDC", `adult lamina epithelial/marginal glial cell`="#9E9AC8", `adult reticular neuropil associated glial cell`="#756BB1", `optic lobe associated cortex glial cell`="#54278F", `cholinergic neuron`="#E41A1C", `gabaergic neuron`="#377EB8", `glutamatergic neuron`="#4DAF4A", `dopaminergic neuron`="#984EA3", `serotonergic neuron`="#FF7F00", `OPN neuron`="#A6CEE3", `columnar neuron T1`="#DEEBF7", `T neuron T3`="#9ECAE1", `T neuron T4/T5`="#3182BD", `distal medullary amacrine neuron Dm3`="#c6e6ba", `distal medullary amacrine neuron Dm8`="#A1D99B", `distal medullary amacrine neuron Dm9`="#31A354", `proximal medullary amacrine neuron Pm4`="#006D2C", `centrifugal neuron C2`="#7FC97F", `centrifugal neuron C3`="#BEAED4", `medullary intrinsic neuron Mi1`="#d8d5eb", `medullary intrinsic neuron Mi4`="#BCBDDC", `medullary intrinsic neuron Mi15`="#756BB1", `transmedullary neuron Tm1`="#b4dee9", `transmedullary neuron Tm2`="#bfdfe7", `transmedullary neuron Tm9`="#7eb5c4", `transmedullary neuron Tm20`="#9fcad5", `transmedullary Y neuron TmY14`="#5da1b3", `transmedullary Y neuron TmY5a`="#3c8ca2", `transmedullary Y neuron TmY8`="#1a7790", `tm5ab`="#014636", `lamina monopolar neuron L1 + L2`="#FDD3A0", `lamina monopolar neuron L3`="#FDBE85", `lamina monopolar neuron L4`="#FD8D3C", `lamina monopolar neuron L5`="#E6550D", `lamina wide-field 2 neuron`="#A63603", `lobular columnar neuron LC12`="#BC80BD", `olfactory receptor neuron`="#CCEBC5", `Poxn neuron`="#FFED6F", `hemocyte`="#999999", `kenyon cell`="#B3E2CD", `gamma kenyon cell`="#FDCDAC", `pigment cell`="#EE82EE")
clusterPalette <- list(`0`="#76C15D", `1`="#539DA4", `2`="#B15928", `3`="#78B69A", `4`="#46A93A", `5`="#EA4041", `6`="#E42123", `7`="#A5D880", `8`="#FDB35B", `9`="#F29A94", `10`="#449F35", `11`="#F98F8F", `12`="#2F83AF", `13`="#A6CEE3", `14`="#84B8D7", `15`="#85C969", `16`="#DBA08D", `17`="#C09B79", `18`="#F68724", `19`="#9CCF90", `20`="#2078B4", `21`="#5D9E43", `22`="#E7392B", `23`="#6D4199", `24`="#E73132", `25`="#E42521", `26`="#F89F5F", `27`="#37A22F", `28`="#8AC395", `29`="#D99B86", `30`="#63A3CB", `31`="#9774B6", `32`="#418EC0", `33`="#CEADC2", `34`="#4190AA", `35`="#F6807F", `36`="#A79C6B", `37`="#F1764A", `38`="#5298C5", `39`="#95D074", `40`="#F37070", `41`="#A48999", `42`="#8F9D5E", `43`="#EEE999", `44`="#AEDC8A", `45`="#FE8B16", `46`="#66A99F", `47`="#EF8D3E", `48`="#56B146", `49`="#FE9324", `50`="#BA9FCC", `51`="#FE8408", `52`="#FEFD97", `53`="#F06060", `54`="#E89458", `55`="#EE6240", `56`="#EB4D36", `57`="#FE9B31", `58`="#769D50", `59`="#7F5999", `60`="#66B952", `61`="#FC8109", `62`="#95C3DD", `63`="#73AED1", `64`="#DBD199", `65`="#C6AED3", `66`="#D5A7A8", `67`="#E29A73", `68`="#8B65AE", `69`="#3183BA", `70`="#AF91C5", `71`="#FBB369", `72`="#FDA33F", `73`="#C9B999", `74`="#FDAB4D", `75`="#7F57A7", `76`="#A382BD", `77`="#927199", `78`="#F48B54", `79`="#E1BF6D", `80`="#FDBB68", `81`="#ED5051", `82`="#F4E889", `83`="#D7AB5F", `84`="#C48243", `85`="#7348A0", `86`="#EAD47B", `87`="#B7A199", `88`="#CD9651", `89`="#BA6D35")
annotatedClusterPalette <- list(`0: unannotated`="#76C15D", `1: unannotated`="#539DA4", `2: unannotated`="#B15928", `3: glutamatergic neuron`="#78B69A", `4: unannotated`="#46A93A", `5: unannotated`="#EA4041", `6: epithelial cell`="#E42123", `7: transmedullary Y neuron TmY8`="#A5D880", `8: cone cell`="#FDB35B", `9: photoreceptor cell`="#F29A94", `10: photoreceptor cell`="#449F35", `11: cone cell`="#F98F8F", `12: cholinergic neuron`="#2F83AF", `13: skeletal muscle of head`="#A6CEE3", `14: gabaergic neuron`="#84B8D7", `15: unannotated`="#85C969", `16: medullary intrinsic neuron Mi4`="#DBA08D", `17: unannotated`="#C09B79", `18: outer photoreceptor cell`="#F68724", `19: epithelial cell`="#9CCF90", `20: skeletal muscle of head`="#2078B4", `21: gabaergic neuron`="#5D9E43", `22: unannotated`="#E7392B", `23: unannotated`="#6D4199", `24: cone cell`="#E73132", `25: T neuron T4/T5`="#E42521", `26: OPN neuron`="#F89F5F", `27: kenyon cell`="#37A22F", `28: outer photoreceptor cell`="#8AC395", `29: olfactory receptor neuron`="#D99B86", `30: unannotated`="#63A3CB", `31: unannotated`="#9774B6", `32: gamma kenyon cell`="#418EC0", `33: unannotated`="#CEADC2", `34: distal medullary amacrine neuron Dm3`="#4190AA", `35: serotonergic neuron`="#F6807F", `36: lamina monopolar neuron L1 + L2`="#A79C6B", `37: adult reticular neuropil associated glial cell`="#F1764A", `38: epithelial cell`="#5298C5", `39: adult optic chiasma glial cell`="#95D074", `40: unannotated`="#F37070", `41: epithelial cell`="#A48999", `42: tm5ab`="#8F9D5E", `43: unannotated`="#EEE999", `44: T neuron T3`="#AEDC8A", `45: unannotated`="#FE8B16", `46: unannotated`="#66A99F", `47: pigment cell`="#EF8D3E", `48: epithelial cell`="#56B146", `49: lobular columnar neuron LC12`="#FE9324", `50: centrifugal neuron C2`="#BA9FCC", `51: Poxn neuron`="#FE8408", `52: transmedullary Y neuron TmY5a`="#FEFD97", `53: lamina monopolar neuron L4`="#F06060", `54: photoreceptor-like cell`="#E89458", `55: lamina wide-field 2 neuron`="#EE6240", `56: columnar neuron T1`="#EB4D36", `57: adult brain cell body glial cell`="#FE9B31", `58: dopaminergic neuron`="#769D50", `59: lamina monopolar neuron L3`="#7F5999", `60: transmedullary neuron Tm1`="#66B952", `61: medullary intrinsic neuron Mi1`="#FC8109", `62: medullary intrinsic neuron Mi4`="#95C3DD", `63: transmedullary neuron Tm9`="#73AED1", `64: adult brain perineurial glial cell`="#DBD199", `65: skeletal muscle of head`="#C6AED3", `66: adult lamina epithelial/marginal glial cell`="#D5A7A8", `67: transmedullary neuron Tm20`="#E29A73", `68: hemocyte`="#8B65AE", `69: adult optic chiasma glial cell`="#3183BA", `70: unannotated`="#AF91C5", `71: transmedullary neuron Tm2`="#FBB369", `72: lamina monopolar neuron L5`="#FDA33F", `73: glutamatergic neuron`="#C9B999", `74: epithelial cell`="#FDAB4D", `75: cholinergic neuron`="#7F57A7", `76: centrifugal neuron C3`="#A382BD", `77: unannotated`="#927199", `78: distal medullary amacrine neuron Dm9`="#F48B54", `79: distal medullary amacrine neuron Dm8`="#E1BF6D", `80: medullary intrinsic neuron Mi15`="#FDBB68", `81: proximal medullary amacrine neuron Pm4`="#ED5051", `82: optic lobe associated cortex glial cell`="#F4E889", `83: unannotated`="#D7AB5F", `84: OPN neuron`="#C48243", `85: unannotated`="#7348A0", `86: unannotated`="#EAD47B", `87: unannotated`="#B7A199", `88: unannotated`="#CD9651", `89: transmedullary Y neuron TmY14`="#BA6D35")

# Random seed
set.seed(42)

Parameters

input_seurat_object <- "./data/Pan_neuro_integrated_FINAL.rds"
output_figure_prefix <- "./figures/Figure_2I/"

I. Reading Seurat object

First step is to read the Seurat objects, previously created

message("Loading Seurat object...")
## Loading Seurat object...
data.seurat <- readRDS(input_seurat_object)
message(ncol(data.seurat), " cells were loaded")
## 23179 cells were loaded
message(nrow(data.seurat), " genes were loaded")
## 23932 genes were loaded

Regulon UMAP plots

p <- FeaturePlot(data.seurat, reduction = "SCENICbinarizedumap", pt.size = 0.5, features = "SCENIC_sima_activating_regulon", order = T, slot = "data") + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_sima_regulon.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_sima_regulon.png"), width = 12, height = 12, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "SCENICbinarizedumap", pt.size = 0.5, features = "SCENIC_sima_activating_regulon_binarized", order = T, slot = "data") + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_sima_regulon_binarized.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_sima_regulon_binarized.png"), width = 12, height = 12, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "SCENICbinarizedumap", pt.size = 0.5, features = "SCENIC_sima_activating_regulon_binarized", order = T, slot = "data") + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_sima_regulon_binarized_small.pdf"), width = 4, height = 4, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_sima_regulon_binarized_small.png"), width = 4, height = 4, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "SCENICbinarizedumap", pt.size = 0.5, features = "SCENIC_crc_activating_regulon", order = T) + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_crc_regulon.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_crc_regulon.png"), width = 12, height = 12, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "SCENICbinarizedumap", pt.size = 0.5, features = "SCENIC_crc_activating_regulon_binarized", order = T) + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_crc_regulon_binarized.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_crc_regulon_binarized.png"), width = 12, height = 12, dpi = 1000, bg = 'white')

Smaller

p <- FeaturePlot(data.seurat, reduction = "SCENICbinarizedumap", pt.size = 0.5, features = "SCENIC_crc_activating_regulon_binarized", order = T) + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_crc_regulon_binarized_small.pdf"), width = 4, height = 4, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "Regulon_UMAP_crc_regulon_binarized_small.png"), width = 4, height = 4, dpi = 1000, bg = 'white')

UMAP plots

p <- FeaturePlot(data.seurat, reduction = "umap_harmony", pt.size = 0.5, features = "SCENIC_sima_activating_regulon", order = T, slot = "data") + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "UMAP_sima_regulon.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "UMAP_sima_regulon.png"), width = 12, height = 12, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "umap_harmony", pt.size = 0.5, features = "SCENIC_sima_activating_regulon_binarized", order = T, slot = "data") + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "UMAP_sima_regulon_binarized.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "UMAP_sima_regulon_binarized.png"), width = 12, height = 12, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "umap_harmony", pt.size = 0.5, features = "SCENIC_sima_activating_regulon_binarized", order = T, slot = "data") + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "UMAP_sima_regulon_binarized_small.pdf"), width = 4, height = 4, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "UMAP_sima_regulon_binarized_small.png"), width = 4, height = 4, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "umap_harmony", pt.size = 0.5, features = "SCENIC_crc_activating_regulon", order = T) + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "UMAP_crc_regulon.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "UMAP_crc_regulon.png"), width = 12, height = 12, dpi = 1000, bg = 'white')
p <- FeaturePlot(data.seurat, reduction = "umap_harmony", pt.size = 0.5, features = "SCENIC_crc_activating_regulon_binarized", order = T) + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "UMAP_crc_regulon_binarized.pdf"), width = 12, height = 12, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "UMAP_crc_regulon_binarized.png"), width = 12, height = 12, dpi = 1000, bg = 'white')

Smaller

p <- FeaturePlot(data.seurat, reduction = "umap_harmony", pt.size = 0.5, features = "SCENIC_crc_activating_regulon_binarized", order = T) + NoLegend()
p

ggsave(p, filename = paste0(output_figure_prefix, "UMAP_crc_regulon_binarized_small.pdf"), width = 4, height = 4, bg = "white")
ggsave(p, filename = paste0(output_figure_prefix, "UMAP_crc_regulon_binarized_small.png"), width = 4, height = 4, dpi = 1000, bg = 'white')